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Abstract 

In  [18],  two  of  the  authors  developed  a  high  order  accurate  numerical  boundary  condition 
procedure  for  hyperbolic  conservation  laws  on  a  Cartesian  mesh,  which  allows  the  compu¬ 
tation  using  high  order  finite  difference  schemes  on  Cartesian  meshes  to  solve  problems  in 
arbitrary  physical  domains  whose  boundaries  do  not  coincide  with  grid  lines.  This  procedure 
is  based  on  the  so-called  inverse  Lax-Wendroff  (ILW)  procedure  for  inflow  boundary  condi¬ 
tions  and  high  order  extrapolation  for  outflow  boundary  conditions.  However,  the  algebra  of 
the  ILW  procedure  is  quite  heavy  for  two  dimensional  (2D)  hyperbolic  systems,  which  makes 
it  difficult  to  implement  the  procedure  for  order  of  accuracy  higher  than  three.  In  this  paper, 
we  first  discuss  a  simplified  and  improved  implementation  for  this  procedure,  which  uses  the 
relatively  complicated  ILW  procedure  only  for  the  evaluation  of  the  first  order  normal  deriva¬ 
tives.  Fifth  order  WENO  type  extrapolation  is  used  for  all  other  derivatives,  regardless  of 
the  direction  of  the  local  characteristics  and  the  smoothness  of  the  solution.  This  makes  the 
implementation  of  a  fifth  order  boundary  treatment  practical  for  2D  systems  with  source 
terms.  For  no-penetration  boundary  condition  of  compressible  inviscid  flows,  a  further  sim¬ 
plification  is  discussed,  in  which  the  evaluation  of  the  tangential  derivatives  involved  in  the 
ILW  procedure  is  avoided.  We  test  our  simplified  and  improved  boundary  treatment  for  Eu¬ 
ler  equations  with  or  without  source  terms  representing  chemical  reactions  in  detonations. 
The  results  demonstrate  the  designed  fifth  order  accuracy,  stability,  and  good  performance 
for  problems  involving  complicated  interactions  between  detonation/shock  waves  and  solid 
boundaries. 
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1  Introduction 


In  [18],  two  of  the  authors  developed  a  high  order  accurate  numerical  boundary  condition 
procedure  for  finite  difference  methods  for  solving  hyperbolic  conservation  laws  on  a  Carte¬ 
sian  mesh,  while  the  physical  domain  can  be  arbitrarily  shaped.  The  challenge  is  that  the 
physical  boundary  does  not  usually  coincide  with  grid  lines  and  can  intersect  the  Cartesian 
grids  in  an  arbitrary  fashion.  For  the  sake  of  accuracy  and  stability,  a  so-called  inverse  Lax- 
Wendroff  (ILW)  procedure  was  proposed.  This  procedure  helps  us  obtain  normal  derivatives 
at  inflow  boundaries  from  time  derivatives  and  tangential  derivatives  by  a  repeated  use  of 
the  partial  differential  equations  (PDEs).  Together  with  high  order  extrapolation  at  outflow 
boundaries,  we  can  impose  accurate  values  for  the  ghost  points  near  the  boundaries  by  a 
Taylor  expansion.  Numerical  examples  in  [18]  show  that  this  method  is  high  order  accurate 
and  stable  under  the  standard  CFL  conditions  determined  by  the  interior  schemes.  In  [19], 
this  approach  was  extended  to  solve  compressible  inviscid  flows  involving  complex  moving 
geometries. 

The  algebra  of  the  ILW  procedure,  which  relies  on  the  PDEs,  can  be  very  heavy  for  fully 
nonlinear  systems  of  equations  in  two  dimensions  (2D).  Of  course,  it  would  be  even  heavier 
for  2D  problems  with  source  terms.  For  this  reason,  in  [18,  19]  we  only  implemented  third 
order  boundary  treatment  for  2D  Euler  equations,  although  our  method  was  designed  to 
achieve  arbitrarily  high  order  of  accuracy  for  general  equations  with  source  terms. 

In  this  paper,  we  aim  to  simplify  and  improve  our  boundary  treatment  procedure  such 
that  it  is  easier  to  implement  for  achieving  higher  order  accuracy.  We  demonstrate  this 
simplified  and  efficient  procedure  on  the  scheme  with  fifth  order  accuracy,  which  is  the  same 
order  as  our  interior  scheme,  for  both  2D  regular  Euler  equations  and  2D  reactive  Euler 
equations  which  contain  source  terms  representing  one-step  chemical  reactions  in  detona¬ 
tions.  We  investigate  each  term  in  the  Taylor  expansion  and  find  that  only  the  first  two 
terms,  i.e.,  the  constant  term  and  the  first  order  normal  derivative  term,  are  crucial  to  be 
implemented  by  the  costly  ILW  procedure  to  ensure  stability.  In  other  words,  near  the  inflow 
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boundaries,  only  the  first  order  normal  derivative  needs  to  be  obtained  by  the  ILW  proce¬ 
dure.  All  the  higher  order  normal  derivatives  can  be  simply  obtained  by  extrapolation  of  a 
suitable  order,  presumably  because  these  terms  are  multiplied  by  0(hk)  for  k  ^  2,  where  h  is 
the  mesh  size,  and  contribute  little  to  stability.  This  simplified  algorithm  significantly  alle¬ 
viates  the  complicated  algebra  in  the  ILW  procedure,  especially  for  2D  systems  with  source 
terms.  We  note  that  a  similar  idea  was  proposed  by  Qiu  and  Shu  for  the  Lax-Wendroff  time 
discretization  [13],  in  which  they  performed  a  weighted  essentially  non-oscillatory  (WENO) 
approximation  only  for  the  reconstruction  of  the  fluxes  to  the  first  order  time  derivative  to 
avoid  spurious  oscillations  and  used  the  inexpensive  central  difference  approximation  for  the 
reconstruction  of  the  higher  order  time  derivatives. 

The  ILW  procedure  can  be  further  simplified  for  no-penetration  boundary  condition  at 
solid  walls  of  compressible  inviscid  flows  with  or  without  chemical  reactions.  In  this  special 
but  important  case,  the  first  time  derivative  of  normal  momentum  can  be  converted  to  the 
first  normal  derivatives  without  introducing  any  tangential  derivatives.  As  a  result,  we  do 
not  need  to  do  any  numerical  differentiation.  The  same  conclusion  can  be  made  if  we  consider 
(reactive)  Euler  equations  in  primitive  variables.  In  fact,  the  effectiveness  and  robustness  of 
imposing  primitive  variables  instead  of  conservative  variables  was  verified  in  [19]  where  the 
no-penetration  boundary  condition  is  prescribed  in  normal  velocity. 

Our  fifth  order  boundary  treatment  requires  fifth  order  extrapolation  for  both  smooth 
solutions  and  solutions  containing  strong  discontinuities.  In  [18],  we  used  Lagrange  extrapo¬ 
lation  or  least  squares  extrapolation  if  the  solution  is  smooth  near  the  boundary.  Otherwise, 
we  developed  a  third  order  accurate  WENO  type  extrapolation  to  prevent  overshoot  or  un¬ 
dershoot  near  the  discontinuities.  We  extend  the  WENO  type  extrapolation  to  fifth  order 
accuracy  in  this  paper  and  use  it  in  both  smooth  and  non-smooth  cases.  In  this  sense,  our 
method  is  improved  to  be  more  uniform. 

There  are  many  successful  boundary  treatment  techniques  based  on  Cartesian  grids  be¬ 
sides  [18].  Unfortunately,  most  of  them  are  at  most  second  order  accurate.  For  example, 
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a  Cartesian  embedded  boundary  method  was  developed  to  solve  the  wave  equation  with 
Dirichlet  or  Neumann  boundary  conditions  in  [8,  9,  10]  and  to  solve  hyperbolic  conservation 
laws  in  [17].  An  internal  boundary  method  based  on  a  reflection  technique  and  interpolation 
was  developed  to  solve  ideal  and  non-ideal  compressible  reactive  flows  in  [25]  and  detona¬ 
tion  problems  in  [11,  21,  22],  Engineering  applications  in  detonations  usually  require  that 
simulations  be  carried  out  for  complex  geometries.  This  involves  interaction  of  detonations 
with  complex  solid  boundaries.  If  the  accuracy  of  numerical  boundary  treatment  is  low,  the 
error  will  pollute  the  computational  results  in  internal  domain  although  the  computational 
accuracy  of  internal  domain  is  of  high  order.  Thus,  it  is  necessary  to  construct  a  high  order 
numerical  boundary  treatment  to  investigate  complicated  interactions  between  detonations 
and  solid  boundaries. 

The  remainder  of  the  paper  is  organized  as  follows.  We  first  give  an  overview  of  the 
discretization  of  the  problem.  Then  we  concentrate  on  the  issue  of  imposing  the  values  of 
ghost  points.  For  ID  Euler  equations,  we  propose  that  the  kth  order  spatial  derivatives,  k  ^ 
2,  can  be  simply  obtained  by  extrapolation  in  all  circumstances.  The  fifth  order  WENO  type 
extrapolation  is  developed  for  this  purpose.  The  improvement  for  ID  problems  is  extended 
to  2D  reactive  Euler  equations.  The  tangential  derivatives  in  the  2D  1LW  procedure  can  be 
avoided  in  the  case  of  no-penetration  boundary  condition  at  solid  walls.  Numerical  examples 
of  (reactive)  Euler  equations  are  shown  in  Section  3.  We  can  see  our  boundary  treatment  is 
fifth  order  accurate,  stable,  and  capable  of  treating  interactions  between  denotation/shock 
waves  and  solid  boundaries  with  good  resolution.  Concluding  remarks  are  given  in  Section 
4. 

2  Scheme  formulation 

We  consider  strongly  hyperbolic  conservation  laws  with  source  term  for  U  =  U(x,y,t)  E  R2 

f  Ut  +  F{U)x  +  G(U)y  =  S{U)  ( x,y)eQ ,  t  >  0, 

\  U(x,y,0)  =  Uo(x,y)  (x,y)eSl, 
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on  a  bounded  domain  ff  with  appropriate  boundary  conditions  prescribed  on  <9f2  at  time  t. 
We  assume  is  covered  by  a  uniform  Cartesian  mesh  with  mesh  size  Ax  =  Ay  =  h,  but  the 
boundary  dtt  does  not  need  to  coincide  with  any  grid  lines.  The  semi-discrete  approximation 
of  (2.1)  is  given  by 

—Ui:j(t)  =  --  (Fi+1/2j  -  Fi_  1/2, ^  (Gi,j+ 1/2  -  ( Uij(t )) ,  (2.2) 

where  Fi+i/2j  and  Gjj+1/2  are  the  numerical  fluxes.  We  use  a  third  order  total  variation 
diminishing  Runge-Kutta  (RK)  method  [16]  to  integrate  the  system  of  ordinary  differential 
equations  (2.2)  in  time 


U 


(1) 

i,j 

(2) 


=  irk  +  Ate  (t/y , 


U)  = 


u 


n+1 

hj 


-un- 

4  i,j 


+  -Uf]  +  -  AtC  ( U(1) 
4  hj  4 


r(2) 


\ul,  +  -3U\:j  +  §Af£  ^ 


(2.3) 


where  £(•)  is  the  operator  defined  by  the  right-hand  side  of  (2.2). 

Suppose  we  have  a  time  dependent  boundary  condition  g(t).  Following  the  idea  in  [3],  we 
can  easily  show  that  for  the  hyperbolic  system  (2.1)  with  source  term  the  following  match 
of  time  levels  at  the  boundary  maintains  the  third  order  accuracy  of  (2.3) 


Un  ~  g(tn), 

C/(1)  ~  g{tn)  +  Atg\tn\  (2.4) 

U (2)  ~  g(tn)  +  tg'{tn)  +  ^AfV(^n)- 

For  simplicity,  we  denote  the  boundary  conditions  for  all  three  stages  at  time  level  t  =  tn  by 
g(tn),  although  g(tn)  is  actually  different  for  each  stage  according  to  (2.4). 

We  use  the  fifth  order  finite  difference  WENO  scheme  with  the  Lax- Friedrichs  flux  split¬ 
ting  [7]  to  form  the  numerical  fluxes  F i+ i/2j-  and  GtJ+l/2  in  (2.2).  The  scheme  requires 
a  seven  point  stencil  in  both  x  and  y  directions.  Near  dVt  where  the  numerical  stencil  is 
partially  outside  of  up  to  three  ghost  points  are  needed  in  each  direction.  We  concentrate 
on  how  to  define  the  values  of  the  ghost  points  at  time  level  t  —  tn  in  the  rest  of  the  paper. 
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2.1  One-dimensional  systems 


We  consider  ID  compressible  Euler  equations  without  source  term 

Ut  +  F(U)X  =  0,  x  G  (—1, 1),  t  >  0, 
where  the  conservative  variables 

‘'-(IMi) 

and  the  flux 

^  U<2  2  \  f  Pu  \ 

F  (U)  —  (7-l)^3  +  Vlq  =  pu2+p  . 

V(Wa  V“(S  +  P)/ 

p,  u,  p  and  E  describe  the  density,  velocity,  pressure  and  total  energy,  respectively.  The 
equation  of  state  has  the  form 

where  7  is  the  (constant)  specific  heat  ratio.  The  sound  speed  is  c  =  \fypfp.  We  consider 
the  right  boundary  x  —  1.  The  left  boundary  can  be  treated  similarly. 

Let  us  discretize  the  interval  (—1,1)  by  a  uniform  mesh 

—1  +  h/2  =  xq  <  x\  <  . . .  <  xn  =  1  —  h/2. 


We  assume  U 0, . . . ,  U n  have  been  updated  from  time  level  tn-\  to  time  level  tn.  Here  we 
suppress  the  tn  dependence  without  causing  any  confusion.  We  proceed  as  in  [18]  to  construct 
values  of  ghost  points  U  n+i,  ■  ■  ■  ,U  tv+3  by  a  fifth  order  Taylor  expansion 

(■ Um)j  =  U™k)’  rn  =  1,2,3,  j  =  N  +  1, . . . ,  N  +  3,  (2.5) 

k= 0 

where  U*rik>  is  a  (5  —  /c)th  order  approximation  of  the  spatial  derivative  yypr 

X=l,  t=tn 

We  first  do  a  local  characteristic  decomposition  to  decide  the  inflow  and  outflow  boundary 


conditions.  Denote  the  Jacobian  matrix  of  the  flux  near  the  boundary  by 


dF(U) 

dU 


u=uN 


A±  (UN) 
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Aj_(U n)  has  three  eigenvalues  Ai  =  ujy  —  cjv,  A2  =  u^,  A 3  =  %  +  cjv  and  a  complete  set  of 
left  eigenvectors  Ii(Un),  h{U n)  which  forms  a  matrix 


L(Un) 


h(UN )  N 

\ 

f  h,l 

h,2 

h,3 

k(uN) 

= 

hi 

h,2 

h,3 

h{UN)  / 

/ 

\  hi 

h,2 

h,3 

For  definiteness,  we  assume  Ai  <  A2  <  0  and  A3  >  0.  Thus  two  boundary  conditions  are 
needed,  for  example,  Um(l,t)  =  gm(t ),  m  =  1,2.  We  define  the  local  characteristic  variables 
Vm  at  grid  points  near  the  boundary  by 


lm{UN)Uj,  m=  1,2,3,  j  =  N 


(2.6) 


We  extrapolate  (Vm)j  to  the  boundary  with  a  fifth  order  WENO  type  extrapolation  and 
denote  the  extrapolated  kth  order  derivative  of  Vm  at  the  boundary  by  V7 nk\  k  =  0, ...  ,4, 
m  =  1,2,3.  The  WENO  type  extrapolation  will  be  described  in  detail  in  the  next  subsection. 

We  impose  U^())  =  gi(tn)  and  U^  0'1  =  g2(tn).  Because  of  the  signs  of  the  corresponding 
eigenvalues,  Vi,  U2  are  the  ingoing  local  characteristic  variables  and  Vi  is  outgoing.  Therefore, 
for  stability  reason,  U3  is  solved  by  the  extrapolation  equation  of  V3 


^3,1^1^  +  h,2U2K'°>  +  ^3,3  Ur 


(0) 


do)  _ 


=  K 


■(o) 


(2.7) 


Next  we  find  the  first  spatial  derivatives  LCP  with  the  ILW  procedure  for  U\  and  f/2,  together 
with  the  extrapolation  equation  of  V3.  Using  the  first  two  equations  of  the  Euler  system,  we 
write  the  first  time  derivatives  of  U\,  U2  as 

dlh  _  _dU2 
dt  dx 

du2  _  dU3  3-7/2  U2dU2  UjdUA 

dt  dx  2  \U\  dx  U\  dx  )  ' 

Notice  that  the  left-hand  side  of  the  above  equations  is  already  known  at  the  boundary.  Um 
can  then  be  solved  by  the  linear  system 
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(2.8) 


(7  -  l)f/*(1)  +  1 


orr*(0) 

2U2  tt*(  1) 
.(0)  ^2 


u;y 


h,l  t/r(1)  +  +  /3,3^ 


f/* (0)  ' 
C/*(0) 
*(1) 


f/; 


*(1) 


T *(1) 


=  K 


*(i) 


Repeatedly  using  the  Euler  equations  and  the  extrapolation  equation,  we  are  able  to  solve 
Umk\  k  =  2,3,4.  However,  the  algebra  of  the  ILW  procedure  becomes  more  and  more 
complicated  for  larger  k.  See  [18]  for  the  formulas  of  k  —  2.  In  fact,  the  scheme  is  still  stable 
if  we  simply  obtain  U*nk\  k  =  2,3,4,  by  extrapolation.  Namely, 


(  ul[k)  \ 


L  (U 


N 


V  u\[k)  ) 


(  v;,k'!  \ 


1 y*(k ) 

V  U*(t)  / 


(2.9) 


In  other  words,  the  complicated  ILW  procedure  is  not  needed  for  spatial  derivatives  of  order 
higher  than  or  equal  to  two,  regardless  of  the  direction  of  the  local  characteristic  variables, 
presumably  because  high  order  derivative  terms  are  multiplied  by  0(hk),  k  =  2,  3, 4,  anyway. 

We  summarize  our  fifth  order  boundary  treatment  at  the  right  boundary  x  =  1  as  follows, 
assuming  that  Uj,  j  =  0, . . . ,  N,  have  been  updated  from  time  level  fn_i  to  time  level  tn. 


1.  Compute  the  eigenvalues  Xm(U n)  and  left  eigenvectors  lm(U n)  of  the  Jacobian  matrix 
A±(U at)  for  m  =  1,2,3.  Decide  the  prescribed  inflow  boundary  conditions  gm(t)  accord¬ 
ing  to  the  signs  of  Am(L/Jv)- 

2.  Form  the  local  characteristic  variables  (Vm)j,  j  =  N  —  4, . . . ,  N,  as  in  (2.6).  Obtain  Vmk\ 
which  is  a  (5  —  k)th  order  approximation  of  the  A; t li  order  derivative  of  Vm  at  the  boundary, 
by  the  WENO  type  extrapolation.  See  the  next  subsection  for  details. 

3.  Solve  for  U*n°\  m  =  1,2,3,  by  the  prescribed  boundary  conditions  and  extrapolated 
values  Vm°\  such  as  (2.7). 


4.  Lise  the  ILW  procedure  to  write  the  first  derivative  of  gm(t)  as  a  linear  combination  of 
first  spatial  derivatives.  Together  with  the  extrapolation  equations,  form  a  linear  system 
with  Um 1 }  as  unknowns,  such  as  (2.8).  Solve  for  Um 1  ^ ,  m  =  1,2,3. 
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*(k) 

5.  For  k  =  2,3,4,  solve  for  Um  by  extrapolation  equations  (2.9).  The  ILW  procedure  is 
not  needed  in  this  step. 

6.  Impose  the  values  of  ghost  points  by  the  Taylor  expansion  (2.5). 

We  finally  emphasize  that  we  use  fifth  order  WENO  type  extrapolation  in  Step  2  regard¬ 
less  of  the  smoothness  of  the  solution. 


2.2  ID  WENO  type  extrapolation 


We  need  to  obtain  a  (5  —  fc)th  order  approximation  of  the  kt\i  order  derivative  of  the  local 
characteristic  variables  Vm  at  the  boundary,  using  the  grid  point  values  in  the  interior  of  the 
domain  (Vm)j,  j  —  N  —  4, . . . ,  N.  See  Step  2  of  the  algorithm  flowchart  in  Section  2.1.  For 
the  ease  of  notation,  we  suppress  the  subscript  of  V  in  this  subsection  since  each  component 
is  extrapolated  by  the  same  method. 

If  V  is  smooth  near  the  boundary,  the  extrapolation  can  be  easily  done  by  constructing 
a  Lagrange  polynomial  p4(x)  interpolating  Vj,  j  =  N  —  4, . . . ,  N,  and  evaluating  the  kth 
order  derivative  of  p4(x)  at  the  boundary.  However,  if  a  discontinuity  appears  in  the  stencil 
S4  =  {xAr_4, . . .  ,xn},  Lagrange  extrapolation  may  lead  to  severe  overshoot  or  undershoot 
near  the  discontinuity.  In  this  situation,  we  prefer  a  lower  order  accurate  but  more  robust 
extrapolation.  A  third  order  accurate  WENO  type  extrapolation  was  developed  in  Section 
2.2  of  [18].  We  extend  it  to  fifth  order  accuracy  here. 

We  consider  the  stencil  S4  as  five  candidate  substencils  Sr  =  {apv-r,  ■  ■  ■ ,  xn},  r  =  0, . . . ,  4. 
On  each  substencil,  we  can  easily  construct  a  Lagrange  polynomial  pr(x)  of  degree  r  such 
that  Pr(xj)  =  Vj,  j  —  N  —  r, . . .  ,N.  Suppose  V  is  smooth  on  S4,  V*^  can  be  extrapolated 
by 


v<k)  = 


dkpr(x ) 


r=0 


dxk 


X=1 


where  d,0  =  h4,  d,\  =  h3,  d2  =  h 2,  d3  =  h,  d4  =  1  —  J0)?=o  dr. 
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We  look  for  WENO  type  extrapolation  in  the  form 


y*(k)  = 


dkpr(x ) 


(2.10) 


where  ur  are  the  nonlinear  weights  depending  on  the  value  of  Vj.  In  the  case  that  V  is 
smooth  in  S4,  we  would  like  to  have 


u0  =  0(h4),  u>i  =  0(h3),  oj2  =  0(h2),  CU3  =  0(h)  and  u4  =  1  — 


(2.11) 


such  that  (2.10)  is  (5  —  /c)th  order  accurate.  The  nonlinear  weights  ur  are  dehned  by 


Ur  = 


Z^S=0  as 


(e  +  (3r)q  ’ 


(2.12) 


where  £  =  10  6,  q  ^  3  is  an  integer  and  [3r  are  the  smoothness  indicators  determined  by 


Po  —  h2 , 

r  rl+h 

Pr  =  Y.  /  ,V 


l-l  (  d 


pr(x )  )  dx,  r  —  1, . . . ,  4. 


We  can  show  that  (2.11)  holds  if  V  is  smooth  in  S4.  If  Sr,  r  ^  1,  contains  a  discontinuity  but 
Sr- 1  does  not,  we  have  us  =  0(hr~l~s ),  s  —  0, . . . ,  r  —  1,  and  us  =  0(h2g+r~1~s),  s  =  r, . . . ,  4. 
Namely,  I/W0-*  reduces  to  a  rth  order  approximation  and  as  h  — >  0  the  weights  assigned  to 
the  non-smooth  substencils  Sr, . . . ,  S4  vanish.  The  proof  is  similar  to  that  in  Section  2.2  of 
[18]  and  is  omitted  here. 

2.3  Two-dimensional  problems 

We  consider  2D  compressible  reactive  Euler  equations 


Ut  +  F(U)x  +  G(U)y  =  S(U),  (x,y)en,  t  >  0, 


(2.13) 
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where 


U  = 


(  Ui\ 

U-2 

U-, 

U4 

\u5J 


G(U)  = 


(  p  \ 

( 

pu 

\ 

pu 

pu-  +  p 

= 

pv 

5 

F(  U)  = 

puv 

E 

u(E+p) 

\PY 

V 

puY  ) 

/ 

pv 

\ 

(  0  ^ 

puv 

0 

pv2  +  p 

,  S(U)  = 

= 

0 

v(E  +  p) 

0 

V 

pvY 

J 

\w  / 

Y  describes  the  reactant  mass  fraction  and  the  source  term  is  assumed  to  be  in  an  Arrhenius 


form 


u  =  -KcpYe~T+/T, 

where  T  =  p/p  is  the  temperature,  T+  is  the  activation  temperature  and  Kc  is  a  constant. 
The  equation  of  state  has  the  form 

E  =  +  L(u2  +  v2)+pQY ,  (2.14) 

7  —  i  z 

where  Q  is  the  (constant)  heat  release  due  to  the  chemical  reaction.  Notice  that  if  Y  =  0, 
then  (2.13)  reduces  to  the  usual  compressible  Euler  equations  without  source  term. 

We  proceed  as  in  Section  2.5  of  [18].  We  assume  the  values  of  all  the  grid  points  inside 
have  been  updated  from  time  level  tn_  1  to  time  level  tn.  For  a  ghost  point  P  =  ( x^y.j ), 
we  find  a  point  Po  =  (xo,yo)  =  on  the  boundary  <911  such  that  the  normal  n(x0)  at  Po 
goes  through  P.  The  sign  of  the  normal  n(x0 )  is  chosen  in  such  a  way  that  it  is  positive  if 
it  points  to  the  exterior  of  fl.  We  set  up  a  local  coordinate  system  at  P0  by 

(•M-r,  rs:)(^)=r(^)’  (2i5) 

where  6  is  the  angle  between  the  normal  n(xo)  and  the  x-axis  and  T  is  a  rotational  matrix. 
The  x-axis  then  points  in  the  same  direction  as  n(x o)  and  the  y-axis  points  in  the  tangential 
direction.  In  this  local  coordinate  system,  the  reactive  Euler  equations  (2.13)  are  written  as 

Ut  +  F(U)£  +  G(U)&  =  S(U),  (2.16) 
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where 


U  = 


/  Pi  \ 

(  p  \ 

U2 

pu 

U3 

= 

pv 

Ua 

E 

\  u5  J 

\PY  ) 

=  T 


For  a  fifth  order  boundary  treatment,  the  value  of  the  ghost  point  P  is  imposed  by  the 
Taylor  expansion 


{U, 


m)i,j 


4  A  k 

V—  /7*(fe) 


fc=0 


///  —  1 . ,5, 


(2.17) 


where  A  is  the  ^-coordinate  of  P  and  UuY  is  a  (5  —  &)th  order  approximation  of  the  normal 

We  assume  U0  is  the  value  of  a  grid  point  nearest  to  Pq  among 


derivative 

9x  \(x,y)=x0,  t=t„ 

all  the  grid  points  inside  hi.  We  denote  the  Jacobian  matrix  of  the  normal  flux  by 


A±(U0)  = 


dF(U) 


dU 


U=U  o 


Aj_(C/0)  has  hve  eigenvalues  Ai  =  uq  —  Co,  A2  =  A3  =  A4  =  ho,  A5  =  ho  +  Co  and  a  complete 
set  of  left  eigenvectors  lm(U0),  rri  —  1, ...  5,  which  forms  a  matrix 


L(U  0)  = 


For  definiteness,  we  assume  Ai  <  0  and  A5  >  A2  =  A3  =  A4  >  0.  Thus  one  boundary  condition 
is  needed  at  P0.  For  example,  the  normal  momentum  is  prescribed  f/2(a3oZ)  =  92 (£)• 

The  local  characteristic  variables  Vm  at  grid  points  near  Pq  are  defined  by 


/  h(U0) 

\ 

(  ^1,1 

h,2 

h,3 

b.4 

h,5 

\ 

h(Uo) 

h,i 

h,2 

h,3 

^2,4 

h,5 

h(Uo) 

= 

h,i 

h,2 

h,3 

h,4 

h,5 

h(Uo) 

U,i 

h,2 

U,3 

U,4 

U,5 

\  h(Uo) 

V  ^5,1 

h,2 

h,3 

h,4 

h,5 

/ 

(Am) /.t ,1/  lm{U 0 ) U ITl  1 ;  •  •  •  j  5,  i±i  Vv)  £  ^i,j i 


(2.18) 


where  C  is  a  set  of  grid  points  used  for  extrapolation.  We  extrapolate  {Vm)^u  to  P0 
and  denote  the  extrapolated  fcth  order  f-derivative  of  Vm  by  Vmk\  k  =  0, . . . , 4.  The  choice 
of  £ij  and  the  fifth  order  WENO  type  extrapolation  in  2D  will  be  described  in  detail  in  the 
next  subsection. 
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We  solve  £/m°')  by  a  linear  system  of  equations 


/  0  1  0  0  0  \ 

/  f>*( 0)  \ 

/  02  (tn )  \ 

^2,1  h,2  h,3  h,4  ^2,5 

U*{0) 

v*(o) 

^3,1  h,2  h,3  h,4  h,5 

f>3(0) 

= 

T/*(0) 

v3 

U,1  h,2  U,3  U,4  U,5 

v*(o) 

\  h,i  h,2  h,3  h,4  h,5  / 

uf>  y 

vU  / 

(2.19) 


Here  the  first  equation  is  the  prescribed  boundary  condition.  The  other  equations  represent 
extrapolation  of  the  four  outgoing  characteristic  variables  Vm,  m  —  2, ...  ,5.  Next,  we  use 
the  1LW  procedure  for  U2.  The  second  equation  of  (2.16)  gives  us 


dU2 

dt 


7 


S£/f +  7_ 


1  ui 


Uf 


U I2 


dUi 

dx 


(3-7) 


+^7  ^fadx  ^  ^  dx+Q ^  ^  dx 


U2dU2 

U i  dx 
dU5  d 


U2U3 


dy  \  Ui 


At  the  boundary,  the  left-hand  side  of  the  above  equation  is  the  known  function  g2(t).  The 
tangential  derivative  on  the  right-hand  side  can  be  computed  by  numerical  differentiation, 
since  we  have  obtained  Unl  '  of  all  the  ghost  points.  Thus  Um  ;  can  be  solved  by  the  linear 
system 


^4*(o) 


/  h;(1)  \ 

fj*W 

u2 

u*{1) 

= 

ul[l) 

uU  / 

-02  (*n)  -  §g 


(0)  fr*(0) 

X3 _ 

W) 


ur’u3 


U- 


V2 

V4 

K 


*(i) 

l(i) 

i 

hi) 

■U) 


\ 


(2.20) 


where 


A*(0)  = 


(  7—3 


U. 


.(0) 


U- 


*(0) 


+  2=1  I 
~  2  I 

U*(07 

/o  \  d*(0) 

(3  7) 

fr*(° ) 

(i  7)~;.<o, 

7-1 

'2,1 

l2,2 

h,3 

^2,4 

' 

'3,1 

h,2 

h,3 

h,4 

4,1 

U,2 

U,3 

U,4 

' 

'5,1 

h,2 

h,3 

h,4 

<5(1-7) 

h,5 

h,5 

l 


l 


4.5 

5.5 


Repeatedly  using  equations  (2.16)  and  the  extrapolation  equations,  we  are  able  to  solve  U1 


*(/c) 


k  =  2,3,4.  However,  the  algebra  of  the  1LW  procedure  becomes  extremely  complicated  for 
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equations  (2.16),  which  makes  the  implementation  impractical.  As  for  ID  problems,  we  can 
avoid  the  cumbersome  ILW  and  use  extrapolation  only.  f/m  ,  k  —  2,3,4,  can  be  solved  by 


L{U  o 


/  Ul[k)  \ 

/  v*{k)  \ 

u*(k) 

V*(k) 

u*(k) 

= 

t r*(k) 

ul{k) 

V*(k) 

uf]  / 

Klt>  / 

(2.2!) 


An  important  class  of  boundary  conditions  for  (reactive)  Euler  equations  is  the  no¬ 
penetration  boundary  condition  at  solid  walls,  i.e.,  u  =  0  or  U2  =  0.  In  this  case,  the 
eigenvalues  Ai  ~  — c0  <  0,  A5  ~  c0  >  0  and  A2  =  A3  =  A4  ~  0.  Since  only  one  boundary 
condition  is  prescribed,  we  consider  Vm,  m  —  2, ...  5,  to  be  outgoing  and  V\  to  be  ingoing, 
which  falls  into  the  same  case  as  discussed  above.  (2.19)  gives  us  =  0.  Then  the  first 
equation  of  (2.20)  reduces  to 


7-1  U\ 


to)' 


U‘tm  +  (1  -  7 (SX"  +  (7  -  W?'  -  <3(7  -  w 


XT, 


‘(0) 


r*(i) 


*(i) 


U\ 


(o) 


r*!1) 


Ui 


(0) 


U', 


■(o) 


RU{ 


(o) 


(2.22) 


where  R  is  the  radius  of  curvature  of  d Vt  at  Pq.  Notice  that  there  is  no  tangential  derivative 
in  (2.22).  Therefore,  we  do  not  need  to  do  any  numerical  differentiation,  which  further 
simplifies  the  implementation. 

(2.22)  can  also  be  derived  by  considering  primitive  variables,  i.e.,  p,u,v,p  and  Y,  in  the 
ILW  procedure.  At  the  solid  walls,  we  obtain 

',2 


dp  l  r 

dx  ^  R  1 


(2.23) 


which  is  equivalent  to  (2.22)  because  of  the  equation  of  state  (2.14).  In  fact,  it  is  some¬ 
times  more  convenient  to  use  primitive  variables  than  to  use  conservative  variables  U.  The 
performance  of  the  former  is  good  as  illustrated  in  [19]. 

We  now  summarize  our  fifth  order  boundary  treatment  for  the  2D  problem  (2.13).  We 
assume  the  values  of  all  the  grid  points  inside  D  have  been  updated  from  time  level  tn_x  to 
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time  level  tn.  Our  goal  is  to  impose  the  value  of  ( Um)i,j ,  m  =  1, . . . ,  5,  for  each  ghost  point 

1.  For  each  ghost  point  ( Xi,yj ),  we  do  the  following  three  steps: 

•  Decide  the  local  coordinate  system  (2.15).  Compute  the  eigenvalues  \m(Uo)  and 
left  eigenvectors  lm(U0)  of  the  Jacobian  matrix  Aj_(U0 )  for  m  —  1, ...  5.  Decide  the 
prescribed  inflow  boundary  conditions  gm(t)  according  to  the  signs  of 

•  Form  the  local  characteristic  variables  (Vm)^ u,  G  £ltJ1  as  in  (2.18).  Extrap- 

*(k) 

olate  [Vm)^u  to  the  boundary  to  obtain  Vm  ,  k  —  0, ...  ,4,  with  fifth  order  WENO 
type  extrapolation.  Details  of  the  2D  extrapolation  will  be  discussed  in  the  next 
subsection. 

•  Solve  for  Um°\  m  —  1, ...  ,5,  by  the  prescribed  boundary  conditions  and  extrapo¬ 
lated  values  Vm°\  such  as  (2.19). 

2.  For  each  ghost  point  ( Xi ,  yj),  use  the  ILW  procedure  to  write  the  first  derivative  of  gm(t) 
as  a  linear  combination  of  first  normal  derivatives  plus  tangential  derivatives.  Together 
with  the  extrapolation  equations,  form  a  linear  system  with  Uffl  as  unknowns,  such  as 
(2.20).  Solve  for  Uml\  m  =  1,  ...5.  For  k  =  2,3,4,  solve  for  U*nk>  by  extrapolation 
equations  (2.21)  where  the  ILW  procedure  is  not  needed. 

3.  Impose  the  values  of  the  ghost  points  by  the  Taylor  expansion  (2.17). 

If  no-penetration  boundary  condition  is  considered  at  solid  walls,  then  the  first  equation 
of  (2.20)  is  replaced  by  (2.22)  in  Step  2  with  other  steps  unchanged.  An  alternative  is  to  use 
primitive  variables  instead  of  conservative  variables. 

2.4  2D  WENO  type  extrapolation 

We  return  to  the  issue  of  2D  fifth  order  WENO  type  extrapolation.  It  is  needed  in  the  second 
bullet  of  Step  1  of  the  algorithm  flowchart  in  Section  2.3,  regardless  of  the  smoothness  of  the 
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solution.  For  stability  reason,  we  used  third  order  least  squares  extrapolation  if  the  solution 
is  smooth  near  the  boundary  or  third  order  WENO  type  extrapolation  otherwise  in  [18]. 
Our  fifth  order  extrapolation  should  combine  the  features  of  both.  Namely,  more  points  are 
needed  to  construct  a  2D  extrapolating  polynomial  and  the  extrapolation  should  not  have 
severe  overshoot  or  undershoot  near  the  discontinuities. 

For  a  given  ghost  point  (xi,  i/j),  we  aim  to  first  construct  a  stencil  £  C  D  for  extrapolation 
and  then  to  obtain  a  (5  —  fc)th  order  approximation  of  ,  which  is  denoted  by 

E*(fe),  k  =  0,...,4.  Compared  with  the  notations  used  in  Section  2.3,  here  we  suppress 
the  subscripts  indicating  the  mth  component  of  V  and  that  the  stencil  £  is  for  the  ghost 
point  ( Xi,yj ).  £  =  (Jr=0  £r  consists  of  five  substencils  £r,  r  —  0, . .  .,4,  each  of  which  is  for 
constructing  a  least  squares  polynomial  pr(x,y )  in  Pr,  i.e., 

Pr(x,y)  =  almXlym 

O^m+l^r 

satisfying 

Pr{x^yu)  =  V^v,  for  all  (xA,2/„)  G  £r.  (2.24) 

Notice  that  pr(x,y )  has  (r  +  l)(r  +  2)/2  degrees  of  freedom.  We  choose  £r  such  that  it 
contains  (r  +  l)2  points  and  (2.24)  holds  in  the  sense  of  least  squares  if  r  >  0. 

We  take  r  —  3  as  an  example  to  show  how  to  choose  £r.  £3  is  composed  of  16  square 
points  sketched  in  Fig.  2.1.  It  contains  four  ID  substencils  l  =  0, . . . ,  3.  Suppose  the 
normal  n(x 0)  (or  x-axis)  intersects  the  grid  line  y  =  yn-i  at  a  point  P*.  We  denote  the  grid 
point  on  the  horizontal  line  y  =  yn-i  which  is  just  to  the  left  of  P;*  by  (xlm,yn-i).  Then  we 
set 

£>1  =  {(xm—i,yn—i),{xm,yn—i),(xrn+i,yn—i),(xm_i_2iyn—i)}r  1  =  0,..., 3, 

and  £3  =  U I0S1.  Notice  that  we  have  to  shift  S0  to  the  left  by  one  grid  point  so  that  S0 
lies  in  £l. 
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Figure  2.1:  The  choice  of  £3  (square  points)  for  2D  WENO  type  extrapolation. 


A)  =  2  K\ 

Pr  =  I  \K\lahl  (DaPr(x,y)f  dxdy,  r  =  1,  —  ,4, 

l^|a|^r  K 

where  a  is  a  multi-index  and  K  =  [x0  —  h/2,  x0  +  h/2]  x  [yQ  —  h/2,  y0  +  h/2].  We  can  show 
that  (2.25)  is  (5  —  &)th  order  accurate  if  V  is  smooth  in  £  and  c or  vanishes  as  h  — *  0  if  £r 
contains  a  discontinuity. 
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3  Numerical  examples 


We  test  our  fifth  order  boundary  treatment  in  this  section.  For  accuracy  tests,  we  take 
At  =  0(h 5/3)  in  the  third  order  RK  time  integration  (2.3)  in  order  to  have  fifth  order 
accuracy  in  time.  For  other  cases,  we  take  At  =  0(h)  and  the  CFL  number  as  0.6.  The 
specific  heat  ratio  7  in  (reactive)  Euler  equations  is  taken  as  7  =  1.4  and  the  power  q  in 
(2.12)  or  (2.26)  is  taken  as  q  =  3,  unless  otherwise  indicated. 

3.1  ID  Euler  equations 

Example  1  We  redo  Example  3(a)  in  [18]  with  the  simplified  and  improved  algorithm.  The 
domain  is  (—77  7r)  and  the  initial  condition  is 

p(x,Q)  =  1  +  0.2  sin  a;, 
u(x,  0)  =  1, 
p(x,  0)  =  2. 

We  want  to  impose  the  boundary  conditions  such  that  the  exact  solution  is  simply  a  trans¬ 
lation  of  the  initial  condition 


p(x,t )  =  1  +  0.2  sin  (x  —  t), 

u(x,t )  =  1, 

p(x,t)  =  2. 

At  both  boundaries,  we  have  Ai  <  0  and  A3  >  A2  >  0.  We  prescribe  p(—7r,t),  u(—7T,t)  and 

p( 7r,  t).  The  density  errors  are  listed  in  Table  3.1.  We  can  clearly  see  the  designed  fifth  order 
convergence. 


Example  2  We  consider  the  interaction  of  two  blast  waves  [24],  The  initial  data  are 


U(x,  0) 


U l  0  <  x  <  0.1, 
U m  0.1  <  x  <  0.9, 
U a  0.9  <  x  <  1, 
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Table  3.1:  Density  errors  and  convergence  rates  of  Example  1.  h  =  27 r/N.  t  =  2. 


N 

L 1  error 

Order 

L°°  error 

Order 

40 

7.02E-05 

5.84E-04 

80 

2.82E-07 

7.96 

2.84E-06 

7.69 

160 

3.16E-09 

6.48 

1.78E-08 

7.31 

320 

9.89E-11 

5.00 

5.58E-10 

5.00 

640 

3.06E-12 

5.02 

1.73E-11 

5.01 

where  pL  =  Pm  =  Pr  =  1 ,  uL  =  uM  =  uR  =  0,  pL  =  103,  pM  =  10  2,  Pr  =  102.  There  are 
solid  walls  at  both  x  =  0  and  x  —  1.  This  problem  involves  multiple  reflections  of  shocks 
and  rarefactions  off  the  walls.  The  density  profile  at  t  —  0.038  is  shown  in  Fig.  3.1(a)  with 
h  =  1/800  and  in  Fig.  3.1(b)  with  h  =  1/1600.  The  reference  solution  is  computed  by 
the  fifth  order  WENO  scheme  with  h  =  1/16,000,  together  with  the  reflection  technique 
at  both  boundaries.  Our  fifth  order  boundary  treatment  gives  an  excellent  non-oscillatory 
resolution.  The  results  are  similar  to  those  obtained  by  a  third  order  boundary  treatment 
in  Example  4  of  [18]. 


(a)  h=  1/800  (b)h  =  1/1600 

Figure  3.1:  The  density  profiles  of  the  blast  wave  problem.  Solid  lines:  reference  solution 
computed  by  the  fifth  order  WENO  scheme  with  h  =  1/16,000,  together  with  the  reflec¬ 
tion  technique  at  boundaries;  Symbols:  numerical  solutions  by  our  fifth  order  boundary 
treatment. 
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3.2  2D  Euler  equations 

Example  3  We  test  the  vortex  evolution  problem  on  a  disk  hi  =  {( x,y )  :  x2  +  y1  <  0.5}. 
The  mean  flow  is  p  —  p—  u  —  v  —  1.  We  add  to  this  mean  flow  an  isentropic  vortex 
perturbation  centered  at  (x0,y0)  =  (0.3,  0.3)  in  (u,v)  and  in  the  temperature  T  =  p/p,  no 
perturbation  in  the  entropy  S  =  p/p1 

(i 5u ,  8v) 

ST 

SS 

where  (x,  y)  =  (x  —  xo,  y  —  y o),  r2  =  x2  +  y2  and  the  vortex  strength  is  e  =  5.  We  regard 
the  exact  solution  U(x,y,t)  as  the  passive  convection  of  the  vortex  with  the  mean  velocity 
and  take  boundary  conditions  from  U(x,y,t)  whenever  needed.  The  number  of  boundary 
conditions  is  determined  by  the  signs  of  the  four  eigenvalues  Xm  which  vary  both  in  space 
and  in  time.  This  problem  was  tested  by  a  third  order  boundary  treatment  in  Example  7  of 
[18].  Here  we  are  able  to  test  our  fifth  order  method,  thanks  to  the  simplified  algorithm.  The 
density  errors  are  listed  in  Table  3.2.  We  can  see  that  the  designed  fifth  order  convergence 
is  achieved  in  the  L 1  norm. 

Table  3.2:  Density  errors  of  the  vortex  evolution  problem  on  a  disk,  t  —  0.1. 


h 

L 1  error 

Order 

L°°  error 

Order 

1/20 

1.94E-07 

2.55E-05 

1/40 

3.66E-09 

5.73 

1.94E-07 

7.04 

1/80 

1.22E-10 

4.91 

1.07E-08 

4.17 

1/160 

4.26E-12 

4.83 

6.00E-10 

4.16 

<7  —  .(1-r3) 


877H 


=  0, 


Example  4  We  test  the  double  Mach  reflection  problem  [24]  which  involves  a  solid  wall  not 
aligned  with  the  grid  lines.  This  problem  is  initialized  by  sending  a  horizontally  moving  shock 
into  a  solid  wall  inclined  by  a  30°  angle.  In  order  to  impose  the  no-penetration  boundary 
condition  by  the  reflection  technique,  people  usually  solve  an  equivalent  problem  that  puts 
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the  wall  horizontal  and  puts  the  shock  60°  angle  inclined  to  the  wall,  see  for  example  [7,  15]. 
We  consider  the  solution  of  the  equivalent  problem  as  our  reference  solution. 

We  solve  the  original  problem  on  a  uniform  mesh  with  our  fifth  order  boundary  treatment 
in  conservative  variables.  The  problem  setup  is  the  same  as  in  Example  8  of  [18].  The 
computational  domain  is  sketched  in  Fig.  3.2(a),  together  with  some  of  the  grid  points  near 
the  wall  indicating  that  the  wall  is  not  aligned  with  the  grid  lines.  Initially  a  right-moving 
Mach  10  shock  is  positioned  at  (0,0).  At  y  —  0,  the  exact  postshock  condition  is  imposed. 
At  the  top  boundary,  the  flow  values  are  set  to  describe  the  exact  motion  of  the  Mach  10 
shock.  Supersonic  inflow  and  outflow  boundary  conditions  are  used  at  the  left  and  right 
boundary,  respectively.  The  power  q  in  (2.26)  is  taken  as  q  =  20  to  stabilize  our  boundary 
treatment.  Fig.  3.2(b)  shows  the  density  contour  with  h  =  1/320  at  t  —  0.2.  A  zoomed- 
in  region  near  the  double  Mach  stem  is  shown  in  Fig.  3.3(a).  The  region  is  rotated  and 
translated  for  the  ease  of  comparison.  In  Fig.  3.3(b),  we  show  the  reference  solution  on  a 
mesh  with  comparable  size.  Figs.  3.3(c)  and  3.3(d)  show  the  density  contours  on  a  refined 
mesh.  Compared  with  the  results  in  [18],  our  results  here  seem  to  be  closer  to  the  reference 
solution  in  terms  of  the  location  of  the  “blown-up”  portion  around  the  double  Mach  stem. 

Example  5  This  is  an  example  involving  a  curved  wall.  The  problem  is  initialized  by  a 
Mach  3  flow  moving  toward  a  circular  cylinder  of  unit  radius  positioned  at  the  origin  on  a  x- 
y  plane.  The  problem  setup  is  the  same  as  in  Example  9  of  [18].  The  computational  domain 
is  the  upper  half  of  the  physical  domain  shown  in  Fig.  3.4(a),  due  to  the  symmetry  of  this 
problem.  At  y  =  0,  the  reflection  technique  is  used.  Supersonic  inflow  boundary  condition  is 
used  at  the  left  boundary  x  =  —  3  and  supersonic  outflow  boundary  condition  is  used  at  the 
top  boundary  y  —  6  and  at  the  right  boundary  x  —  0.  Our  fifth  order  boundary  treatment  in 
primitive  variables  with  q  =  10  in  (2.26)  is  applied  at  the  surface  of  the  cylinder.  Notice  that 
although  the  initial  condition  is  incompatible  with  the  no-penetration  boundary  condition, 
we  do  not  encounter  any  problem  if  using  primitive  variables,  while  difficulty  arises  if  we  use 
conservative  variables  as  reported  in  [18]. 
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2.5- 


Figure  3.2:  Left:  The  computational  domain  (solid  line)  of  the  double  Mach  reflection 
problem.  The  dashed  line  indicates  the  computational  domain  used  in  [7,  15].  The  square 
points  indicate  some  of  the  grid  points  near  the  wall.  Illustrative  graph,  not  to  scale.  Right: 
Density  contour  of  double  Mach  reflection,  30  contour  lines  from  1.731  to  20.92.  h  =  1/320. 


(a)  h  =  1/320,  original  problem  (b)  h  =  -s/3/480,  equivalent  problem 


Figure  3.3:  Density  contours  of  double  Mach  reflection,  30  contour  lines  from  1.731  to 
20.92.  Zoomed-in  near  the  double  Mach  stem.  The  plots  in  the  left  column  are  rotated  and 
translated  for  comparison. 
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The  pressure  contour  at  steady  state  is  shown  in  Figs.  3.4(b)  and  3.4(c)  with  different 
mesh  sizes.  The  bow  shock  is  well-captured  by  our  method.  For  a  more  quantitative  verifica¬ 
tion,  we  take  advantage  of  the  entropy  along  the  surface,  which  can  be  computed  analytically 
by  using  the  Rankine-Hugoniot  conditions  for  the  streamline  normal  to  the  bow  shock  at 
y  —  0.  Since  there  is  usually  no  grid  point  located  on  the  surface,  we  compute  the  entropy 
errors  of  the  state  Um°\  rn  =  1, . .  .4,  which  is  the  constant  term  in  the  Taylor  expansion 
(2.17),  for  all  the  ghost  points.  We  can  see  superlinear  convergence  rates  in  Table  3.3. 


(b)  k  =  1/40 


(c)  h  =  1/80 


Figure  3.4:  Left:  Physical  domain  of  flow  past  a  cylinder.  The  computational  domain  is  the 
upper  half.  The  square  points  indicate  some  of  the  grid  points  near  the  cylinder.  Illustrative 
sketch,  not  to  scale.  Others:  Pressure  contours  of  flow  past  a  cylinder,  20  contour  lines  from 
2  to  15. 


Table  3.3:  L°°  entropy  errors  on  the  surface  of  the  cylinder  and  rates  of  convergence  of 
Example  5. 


h 

L°o  error 

Order 

1/20 

1.68E-02 

1/40 

7.06E-03 

1.25 

1/80 

1.88E-03 

1.91 

1/160 

7.50E-04 

1.33 
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3.3  2D  reactive  Euler  equations 


Example  6  We  would  like  to  test  the  accuracy  of  our  boundary  treatment  for  2D  reactive 
Euler  equations  in  this  example.  Unfortunately,  it  is  difficulty  to  construct  an  analytical 
solution  with  no-penetration  boundary  condition  at  solid  curved  walls.  We  follow  the  idea 
in  Example  4  of  [19]  to  construct  isentropic  flows  so  that  we  are  able  to  measure  the  entropy 
errors.  We  set  the  source  term  uj  —  0  in  (2.13)  such  that  isentropic  flows  can  be  maintained 
as  long  as  the  flows  keep  smooth.  Namely,  we  only  consider  the  transport  of  the  chemical 
reactants.  The  heat  release  is  set  to  Q  =  2. 

We  consider  a  region  [—4, 4]  x  [—4, 4]  with  all  the  boundaries  as  solid  walls.  A  circular 
cylinder  with  unit  radius  is  centered  at  the  origin.  The  computational  domain  is  then 
=  [—4,4]  x  [—4,4]  D  {(x,y)  :  x2  +  y2  >  1}.  The  initial  conditions  should  be  compatible 
with  the  no-penetration  boundary  condition  at  the  solid  walls  and  at  the  surface  of  the 
cylinder.  For  this  purpose,  we  set 


where 


u(x,y,  0) 
v{x,y,  0) 


Ai  {x,y)ui{x,y), 

\  M%,y)vi(x,y)  +  A 2(x,y)v2(x,y)] , 


Ai  (x,y) 

Mxiy) 

ui  (x,y) 
vi(x,y) 
v2(x,y) 


(4a/2-1) 

^^16  +  y2  —  lj  (a/16  +  x 2  —  l) 
\J  x2  +  y2{x2  —  16)  (y2  —  16) 

(^4?  -  ie)  (A?  - 16) ' 

sill  (  X.  j  sill2  (Xjj  , 

ui(y,x), 

x2  +  y2-l )  sin(^)  . 


The  other  initial  conditions  are  p(x,  y,  0)  =  p(x,  y,  0)  = 

1 


Y(x,y,  0)  = 


2  L 


1  +  sinz  (  —x 


1  and 
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We  apply  our  fifth  order  boundary  treatment  in  primitive  variables  at  the  surface  of  the 
cylinder  and  the  reflection  technique  at  the  solid  walls.  Fig.  3.5  shows  the  density  contour 
with  h  =  1/20  at  t  —  0.5.  The  entropy  errors  at  the  same  time  level  are  listed  in  Table  3.4. 
Fifth  order  convergence  is  achieved. 

4 
3 
2 
1 

>  0 

-1 
-2 
-3 
-4 

Figure  3.5:  Density  contour  of  Example  6.  h  =  1/20,  t  =  0.5. 


Table  3.4:  Entropy  errors  and  convergence  rates  of  Example  6.  t  =  0.5. 


h 

L 1  error 

Order 

L°°  error 

Order 

1/5 

1.44E-03 

5.67E-04 

1/10 

5.00E-05 

4.84 

2.72E-05 

4.38 

1/20 

8.66E-07 

5.85 

5.03E-07 

5.76 

1/40 

2.26E-08 

5.26 

1.81E-08 

4.80 

Example  7  We  start  to  consider  problems  involving  detonations.  It  is  a  rapid  regime  of 
burning  in  which  a  strong  shock  ignites  the  combustible  mixture  of  gases  and  the  burning 
proceeds  to  equilibrium  behind  the  shock,  while  the  energy  released  continues  to  help  to 
drive  the  shock  [14] .  The  first  problem  is  the  reflection  of  detonation  waves  on  a  wedge,  see 
[6,  20]  for  experiments  and  [12,  23]  for  numerical  simulations.  The  computational  domain 
is  sketched  in  Fig.  3.6(a).  A  Zeldovich,  von  Neumann  and  Doering  (ZND)  profile  with 
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T+  =  50,  Q  =  50,  /  =  1.0,  7  =  1.2  is  initially  located  at  x  =  —118.2.  Here  /  is  the  overdrive 
factor  defined  by  /  =  D2  /  Dq j,  where  D  is  the  detonation  velocity  and  Hcj  is  the  velocity 
of  the  corresponding  Chapman- Jouguet  wave.  The  constant  prestate  is  ( p,u,v,p,Y )  = 
(1,0,0, 1,1).  See  [2,  26]  for  details  about  the  ZND  profile.  The  density  and  pressure  fields 
of  this  profile  are  perturbed  so  that  the  instability  of  the  flow  results  in  the  development  of 
triple  points,  whose  trajectories  form  a  cellular  pattern  [2,  4,  14],  The  perturbation  is  in  the 
region  [—117.0,  —115.8]  x  [0, 100]  in  which 


p(x,y,  0)  =p(x,y,  0)  = 


1.8  if  y  e  [20,40]  U  [60,; 
0.2  otherwise. 


Supersonic  inflow  and  outflow  boundary  conditions  are  used  at  x  =  —120  and  at  x  =  120, 
respectively.  The  reflection  technique  is  applied  at  y  =  0  and  at  y  =  100.  Our  fifth  order 
boundary  treatment  in  conservative  variables  with  q  =  40  in  (2.26)  is  used  at  the  surface  of 
the  wedge. 


(a)  Example  7  (b)  Example  8 

Figure  3.6:  Left:  Computational  domain  of  Example  7.  Right:  Physical  domain  of  Example 
8.  The  computational  domain  is  the  upper  half.  Both  are  Illustrative  sketches,  not  to  scale. 


We  first  set  the  wedge  angle  6  =  19.3°.  We  run  the  code  until  t  =  28  with  h  =  0.15. 
Fig.  3.7(a)  shows  the  cellular  detonation  pattern  determined  by  recording  the  history  of  the 
maximum  pressure.  The  cellular  pattern  is  similar  to  the  experimental  results  in  Figs.  5, 
10,  11  of  [6]  and  the  numerical  results  in  Fig.  8  of  [23].  From  Fig.  3.7(a),  we  can  observe 
the  main  feature  of  the  detonation  wave  reflection.  There  is  a  sharp  dividing  line,  made  up 
by  a  triple  point  trajectory,  emerging  near  the  tip  of  the  wedge  and  extending  downstream. 
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The  detonation  cells  below  the  line  are  smaller  than  those  above  the  line  and  are  distorted 


in  shape. 


(a)  0  =  19.3°  (b)  6  =  30° 

Figure  3.7:  Cellular  patterns  of  detonation  waves  over  a  wedge  of  angle  6. 


Next,  we  increase  the  wedge  angle  to  6  =  30°  and  keep  other  parameters  unchanged. 
The  cellular  pattern  at  t  =  30  with  h  =  0.15  is  shown  in  Fig.  3.7(b).  As  the  experimental 
results  in  Fig.  1(d)  of  [20],  Fig.  7  of  [6]  and  the  numerical  result  in  Fig.  9  of  [12],  the  cellular 
pattern  is  not  clear  in  this  case  since  the  spatial  scale  of  the  region  between  the  wedge  and 
the  reflection  of  transverse  waves  is  small.  However,  the  dividing  line  is  still  quite  sharp. 

Example  8  We  consider  the  shock  focusing  problem.  The  physical  domain  is  sketched  in 
Fig.  3.6(b).  The  curve  represents  a  solid  wall  with  parabolic  shape.  A  right-moving  Mach 
2.5  shock  is  initially  positioned  at  x  —  —dw,  where  dw  is  the  depth  of  the  parabolic  wall. 
The  preshock  state  is  (p,  u,  v,p,  Y )  =  (0.2,  0,  0,  0.2, 1).  As  the  incident  shock  wave  enters  the 
cavity  of  the  end-wall,  there  are  local  zones  with  elevated  pressure  and  temperature.  The 
formation  of  such  hot  zones  causes  self-ignition  which  is  sometimes  followed  by  detonation 
initiation.  See  [5]  for  experiments  and  [1]  for  numerical  simulations. 

We  set  T+  =  20,  Q  =  50  and  Kc  =  30.  The  computational  domain  is  the  upper  half  of 
the  physical  domain  due  to  the  symmetry  of  this  problem.  The  reflection  technique  is  used 
at  y  =  0  and  at  y  =  2.7.  Supersonic  inflow  boundary  condition  is  used  at  the  left  boundary 
x  =  —5.  Our  fifth  order  boundary  treatment  in  primitive  variables  with  q  =  40  in  (2.26)  is 
applied  at  the  parabolic  wall. 
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In  the  first  case,  we  consider  a  parabola  y 2  =  —2.7a;  with  dw  =  2.7.  The  temperature 
contours  with  h  =  1/160  and  h  =  1/320  at  different  times  are  shown  in  Fig.  3.8.  At  t  =  0.8, 
a  high  temperature  ignition  zone  is  formed  at  the  lateral  surface  near  the  Mach  stem  (see  the 
top  row  of  Fig.  3.8),  which  shows  that  focusing  induces  exothermic  reaction  of  combustible 
gases.  The  primary  ignition  starts  and  the  developing  detonation  waves  are  focused  around 
the  center  line  (see  the  middle  row).  Similar  structure  persists  until  the  detonation  leaves 
the  cavity  of  the  end- wall  (see  the  bottom  row). 

In  the  second  case,  we  consider  a  shorter  parabola  y 2  =  —4.5a;  with  dw  =  1.62.  Other 
parameters  are  kept  the  same.  The  temperature  contours  with  h  =  1/160  and  h  =  1/320 
at  different  times  are  shown  in  Fig.  3.9.  We  again  obtain  convergent  solutions  with  good 
resolutions. 

4  Concluding  remarks 

We  improve  the  high  order  numerical  boundary  condition  developed  in  [18]  in  two  aspects. 
First,  we  find  that  the  ILW  procedure  is  not  needed  for  obtaining  kt\i  order  normal  derivative 
for  k  ^  2,  regardless  of  the  direction  of  the  local  characteristics.  This  makes  the  implementa¬ 
tion  of  fifth  order  boundary  treatment  practical  for  2D  nonlinear  systems  with  source  terms. 
For  no-penetration  boundary  condition  at  solid  walls,  the  tangential  derivative  terms  involved 
in  the  ILW  procedure  can  be  avoided  to  further  simplify  the  implementation.  Secondly,  a 
fifth  order  WENO  type  extrapolation  is  designed  such  that  our  algorithm  keeps  the  same 
for  solutions  with  or  without  discontinuities.  Various  numerical  examples,  with  or  without 
chemical  reactions,  demonstrate  that  our  improved  boundary  treatment  is  fifth  order  accu¬ 
rate  and  performs  well  even  if  there  are  complicated  interactions  between  detonation/shock 
waves  and  solid  boundaries.  In  future  work,  we  will  try  to  extend  this  boundary  treatment  to 
high  order  discontinuous  Galerkin  methods  for  compressible  inviscid  flows  involving  complex 
geometries  on  rectangular  meshes. 


Figure  3.8:  Temperature  contours  of  the  shock  focusing  problem.  29  contour  lines  in  the 
respective  range.  dw  =  2.7.  The  color  outside  the  computational  domain  is  not  relevant. 
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Figure  3.9:  Temperature  contours  of  the  shock  focusing  problem.  29  contour  lines  in  the 
respective  range.  dw  =  1.62.  The  color  outside  the  computational  domain  is  not  relevant. 


30 


Acknowledgments 

The  research  of  S.  Tan  is  supported  by  AFOSR  grant  FA9550-09-1-0126.  The  research  of 
C.  Wang  is  supported  by  New  Century  Excellent  Talents  in  University  under  grant  num¬ 
ber  NCET-08-0043,  NSFC  grant  10972040  and  the  Foundation  of  State  Key  Laboratory  of 
Explosion  Science  and  Technology  (Grant  No.  ZDKT11-01).  The  research  of  C.-W.  Shu  is 
supported  by  AFOSR  grant  FA9550-09-1-0126  and  NSF  grant  DMS-1112700.  The  research 
of  J.  Ning  is  supported  by  NSFC  grant  11032002. 

References 

[1]  A.M.  Bartenev,  S.V.  Khornik,  B.E.  Gelfand,  H.  Gronig  and  H.  Olivier,  Effect  of  re¬ 
flection  type  on  detonation  initiation  at  shock-wave  focusing,  Shock  Waves,  10  (2000), 
205-215. 

[2]  A.  Bourlioux  and  A.J.  Majda,  Theoretical  and  numerical  structure  of  unstable  detona¬ 
tions,  Philosophical  Transactions  of  the  Royal  Society  of  London  A,  350  (1995),  29-68. 

[3]  M.H.  Carpenter,  D.  Gottlieb,  S.  Abarbanel  and  W.-S.  Don,  The  theoretical  accuracy  of 
Runge-Kutta  time  discretizations  for  the  initial  boundary  value  problem:  a  study  of  the 
boundary  error,  SIAM  Journal  on  Scientific  Computing,  16  (1995),  1241-1252. 

[4]  V.N.  Gamezo,  D.  Desbordes  and  E.S.  Oran,  Two-dimensional  reactive  flow  dynamics 
in  cellular  detonation  waves ,  Shock  Waves,  9  (1999),  11-17. 

[5]  B.E.  Gelfand,  S.V.  Khornik,  A.M.  Bartenev,  S.P.  Medvedev,  H.  Gronig  and  H.  Olivier, 
Detonation  and  deflagration  initiation  at  the  focusing  of  shock  waves  in  combustible 
gaseous  mixture ,  Shock  Waves,  10  (2000),  197-204. 

[6]  C.  Guo,  D.  Zhang  and  W.  Xie,  The  Mach  reflection  of  a  detonation  based  on  soot  track 
measurements,  Combustion  and  Flame,  127  (2001),  2051-2058. 


31 


[7]  G.-S.  Jiang  and  C.-W.  Shu,  Efficient  implementation  of  weighted  ENO  schemes ,  Journal 
of  Computational  Physics,  126  (1996),  202-228. 

[8]  H.-O.  Kreiss  and  N.A.  Petersson,  A  second  order  accurate  embedded  boundary  method 
for  the  wave  equation  with  Dirichlet  data,  SIAM  Journal  on  Scientific  Computing,  27 
(2006),  1141-1167. 

[9]  H.-O.  Kreiss,  N.A.  Petersson  and  J.  Ystrom,  Difference  approximations  for  the  second 
order  wave  equation,  SIAM  Journal  on  Numerical  Analysis,  40  (2002),  1940-1967. 

[10]  H.-O.  Kreiss,  N.A.  Petersson  and  J.  Ystrom,  Difference  approximations  of  the  Neumann 
problem  for  the  second  order  wave  equation,  SIAM  Journal  on  Numerical  Analysis,  42 
(2004),  1292-1323. 

[11]  J.-G.  Ning  and  L.-W.  Chen,  Fuzzy  interface  treatment  in  Eulerian  method ,  Science  in 
China  Series  G:  Physics,  Mechanics  &  Astronomy,  47  (2004),  550-568. 

[12]  S.  Ohyagi,  T.  Obara,  F.  Nakata  and  S.  Hoshi,  A  numerical  simulation  of  reflection 
processes  of  a  detonation  wave  on  a  wedge,  Shock  Waves,  10  (2000),  185-190. 

[13]  J.  Qiu  and  C.-W.  Shu,  Finite  difference  WENO  schemes  with  Lax-Wendroff-type  time 
discretizations,  SIAM  Journal  on  Scientific  Computing,  24  (2003),  2185-2198. 

[14]  G.J.  Sharpe  and  S.A.E.G.  Falle,  Two-dimensional  numerical  simulations  of  idealized 
detonations,  Philosophical  Transactions  of  the  Royal  Society  of  London  A,  456  (2000), 
2081-2100. 

[15]  J.  Shi,  Y.-T.  Zhang  and  C.-W.  Shu,  Resolution  of  high  order  WENO  schemes  for  com¬ 
plicated  flow  structures,  Journal  of  Computational  Physics,  186  (2003),  690-696. 

[16]  C.-W.  Shu  and  S.  Osher,  Efficient  implementation  of  essentially  non- os  dilatory  shock¬ 
capturing  schemes,  Journal  of  Computational  Physics,  77  (1988),  439-471. 


32 


[17]  B.  Sjogreen  and  N.A.  Petersson,  A  Cartesian  embedded  boundary  method  for  hyperbolic 
conservation  laws,  Communications  in  Computational  Physics,  2  (2007),  1199-1219. 

[18]  S.  Tan  and  C.-W.  Shu,  Inverse  Lax-  Wendroff  procedure  for  numerical  boundary  condi¬ 
tions  of  conservation  laws,  Journal  of  Computational  Physics,  229  (2010),  8144-8166. 

[19]  S.  Tan  and  C.-W.  Shu,  A  high  order  moving  boundary  treatment  for  compressible  inviscid 
flows,  Journal  of  Computational  Physics,  230  (2011),  6023-6036. 

[20]  G.O.  Thomas  and  R.L1.  Williams,  Detonation  interaction  with  wedges  and  bends,  Shock 
Waves,  11  (2002),  481-492. 

[21]  C.  Wang,  T.-B.  Ma  and  J.  Lu,  Influence  of  obstacle  disturbance  in  a  duct  on  explosion 
characteristics  of  coal  gas,  Science  China:  Physics,  Mechanics  &  Astronomy,  53  (2010), 
269-278. 

[22]  C.  Wang,  T.-B.  Ma  and  J.-G.  Ning,  Tracking  method  for  multi-material  interfaces  and 
its  application  in  shaped  charge ,  Journal  of  Computational  and  Theoretical  Nanoscience, 
5  (2008),  1512-1516. 

[23]  G.  Wang,  D.  Zhang,  K.  Liu  and  J.  Wang,  An  improved  CE/SE  scheme  for  numerical 
simulation  of  gaseous  and  two-phase  detonations,  Computers  &  Fluids,  39  (2010),  168- 
177. 

[24]  P.  Woodward  and  P.  Colclla,  The  numerical  simulation  of  two-dimensional  fluid  flow 
with  strong  shocks,  Journal  of  Computational  Physics,  54  (1984),  115-173. 

[25]  S.  Xu,  T.  Aslam  and  D.S.  Stewart,  High  resolution  numerical  simulation  of  ideal  and 
non-ideal  compressible  reacting  flows  with  embedded  internal  boundaries,  Combustion 
Theory  and  Modelling,  1  (1997),  113-142. 


33 


[26]  Z.-C.  Zhang,  S.T.J.  Yu,  H.  He  and  S.-C.  Chang,  Direct  calculations  of  two-  and  three- 
dimensional  detonations  by  an  extended  CE/SE  method ,  AIAA  Paper  No.  2001-0476, 
2001. 


34 


